Stabilization of BEC droplet in free space by feedback control of interatomic 

interaction 
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A self-trapped Bose-Einstein condensate in three-dimensional free space is shown to be stabi- 
lized by feedback control of the interatomic interaction through nondestructive measurement of the 
condensate's peak column density. The stability is found to be robust against poor resolution and 
experimental errors in the measurement. 
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I. INTRODUCTION 

Matter-wave bright solitons of a Bose-Einstein con- 
densate (BEC) in a quasi-one dimensional (ID) trapping 
potential have been realized by the ENS group [l| and 
the Rice group @] . The stability of this self-trapped state 
is achieved by the balance between the attractive inter- 
atomic interaction and the quantum kinetic pressure. In 
ID, the bright soliton is stable and robust against noise. 
In 2D and 3D, however, the balance between the attrac- 
tive interaction and the kinetic pressure is precarious, 
and infinitesimal deviations from the stationary state 
cause a collapse or expansion of the system. While self- 
trapped liquids, such as a raindrop, are quite common, a 
self-trapped gas is a novel state of matter. Such a stable 
self-trapped state in a gaseous BEC in 2D or 3D, if it 
can be realized, might be referred to as a matter-wave 
droplet or simply a BEC droplet. 

A scheme proposed in Refs. 0, to stabilize a BEC 
droplet is to oscillate the interatomic interaction rapidly 
using Feshbach resonance E| . The rapid oscillation of 
the interaction produces an effective potential that pre- 
vents the condensate from collapsing. This phenomenon 
is similar to the stabilization of an inverted pendulum 
by an oscillating pivot 0. Several researchers have 
demonstrated the stabilization of a BEC droplet in 2D 
with oscillating interactions by numerically solving the 
Gross-Pitaevskii (CP) equation dHHH- Montesinos 
et al. ^3 have shown that this stabilization in 2D is 
also possible for a multicomponent BEC. Matuszewski et 
al. [llj have predicted 3D breather solitons confined in 
a ID lattice. However, in 3D free space, it appears that 
an oscillating interaction alone cannot stabilize a BEC 
droplet due to dynamical instabilities |j, |8j . By taking 
into account the effect of energy dissipation, which always 
exists in realistic situations, we have shown that a BEC 
dro plet with oscillating interactions can be stabilized in 

3D m. 

In the present paper, we show that a BEC droplet in 
3D free space can be stabilized by feedback control of 
the interaction through nondestructive measurement of 
the condensate's peak column density. Real-time moni- 



toring of the density profile of a condensate is possible by 
using a nondestructive in-situ imaging method 13] . The 
collapse or expansion of the condensate can be prevented 
by a decrease or increase, respectively, in the strength of 
the attractive interaction if the peak density of the con- 
densate increases above or decreases below a prescribed 
value. Thus, the shape of the BEC droplet can be main- 
tained without collapse or expansion by negative feed- 
back from the result of the real-time measurement. We 
will also examine the effect of experimental imperfections 
in the real-time measurement, such as spatial and time 
resolutions and experimental errors, and show that the 
stabilized BEC droplet is robust against these imperfec- 
tions. 

This paper is organized as follows. Section UTI discusses 
the stationary solutions of the GP equations in ID, 2D, 
and 3D free space. Section II I II numerically investigates 
the dynamics of the condensate under feedback control, 
and shows that a BEC droplet in 3D can be stabilized 
for a wide range of parameters. Section IIVI presents 
variational analysis using a Gaussian trial function. Sec- 
tion [V] studies the stability against measurement errors 
and finite measurement resolution, and shows that a BEC 
droplet is robust against these experimental imperfec- 
tions. Finally, Sec. I VII concludes this paper. 



II. STATIONARY SOLUTIONS OF THE 
GROSS-PITAEVSKII EQUATION IN FREE 
SPACE 

We first consider a stationary state of a BEC with an 
attractive interaction in free space. The dynamics of the 
BEC are described by the GP equation, 



ih 



~dt 



2m 



v> + Vip - 



Airfi 2 a 



MV, (i) 



where m is the mass of an atom, V is the external po- 
tential, and a is the s-wave scattering length. The wave 
function is normalized with J \ tp\ 2 dr — N, with N being 
the number of atoms. 

When the external potential is given by V = 



id 



(x 2 +y 2 )/2 and Tiwid is much larger than the other 



2 



characteristic energies, the degrees of freedom of the BEC 
in the x and y directions are frozen and the system be- 
haves as an effective ID system. Writing the wave func- 
tion as 



and integrating Eq. over x and y, we obtain an effec- 
tive ID equation: 

di) z h 2 d 2 ^ z 2 
^ = "2^^ +ffld| ^ 1 ^ (3) 

where gid = 2hujida. If a < 0, this equation has a well- 
known soliton solution |14|: 



4>z 



(4) 



where /i = —m,u) 2 d a 2 N 2 /2, rj = muiid\a\N/h, and z$ is 
the position of the peak. This is the ground state so- 
lution, and it has no dynamical instabilities. Thus an 
attractive BEC is stable in ID. 



When V = 



-2d 



z 2 /2 and hu>2d is much larger than 



the other characteristic energies, the system behaves as 
an effective 2D system. Substituting the wave function 



e-^* 2 -^^ xy {x lV ) (5) 



into Eq. Q), and integrating the result over z, we obtain 
an effective 2D equation: 



ih d^y ^ (d 2 | d 2 
dt 2m \dx 2 dy 2 



+ 92d\lpxy\ 2 '4'xv, (6) 



where g2d — {&n'h i a 2 LU2d/'rn) 1 ^ 2 ■ At the critical strength 
of the interaction 
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Eq. © has a stationary self-trapped solution, which is 
known as the Townes soliton [lj|. If the peak of the 
Townes soliton is located at the origin r = 0, the wave 
function is axisymmetric and the density \ijj\ 2 monoton- 
ically decreases to zero for r — > oo. The Townes soliton 
also has a scaling property; if ipxy(x, y) is a stationary so- 
lution of Eq. 10, the scaled wave function aip xy (ax,ay) 
with an arbitrary scaling parameter a is also a stationary 
solution of Eq. ©. However, the Townes soliton is dy- 
namically unstable in the sense that an infinitesimal de- 
viation from the stationary solution grows exponentially 
in time, and therefore the Townes soliton eventually col- 
lapses or expands. 

When V = 0, i.e., in 3D free space, Eq. (Q) with a < 
has a stationary self-trapped solution that is dynamically 
unstable like the Townes soliton. The density of this sta- 
tionary solution has spherical symmetry and monotoni- 
cally decreases to zero for r — > oo ^(|. The striking dif- 
ference between this 3D stationary state and the Townes 
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FIG. 1: Stationary unstable solution of Eq. iJHJ (solid curve) 
for g = — (27r) 3//2 . The Gaussian function e~ r /ty 3 ^ 2 is super- 
imposed as a dashed curve for comparison. See Sec.lTvlfor 
the choice of parameters. The inset shows the same functions 
on a logarithmic scale to show the difference in the large-r 
behavior. 



soliton is in their scaling properties. Normalizing the 
length, time, and wave function in units of £, m£ 2 /h, and 
(N/i 3 ) 1 / 2 , respectively, where I is an arbitrary length 
scale, Eq. reduces to 



V2 / a. 1/12/ 



(8) 



where g = AirNa/i. Therefore, if ip(r,t) is a solution of 
the GP equation with scattering length a, the scaled wave 
function c?l 2 ^>{ar^ a 2 t) also satisfies the GP equation 
with scattering length a /a. This indicates that there 
always exists an unstable stationary solution for any g < 
0, and the solutions for different <?'s are related by the 
scaling transformation. In contrast, in 2D, the Townes 
soliton exists only for a particular value of <?2d = g^d ( see 
Eq. Q) for an arbitrary scaling parameter a. 

In Fig. 2] we plot the density profile of the unstable 
stationary state in 3D for g = — (2-7r) 3 / 2 , which is numeri- 
cally obtained by the Newton-Raphson method 01 • The 
inset of Fig. ^ shows the logarithmic plot of the large-r 
behavior. We find that the tail of the unstable station- 
ary state (solid curve) is longer than that of the Gaus- 
sian wave function e~ r /tt 3 / 2 (dashed curve). Generally, 
in d dimensions the unstable stationary wave function 
has an asymptotic form r ( 1 - d )/ 2 e - cr ' with constant c for 
r — ► oo |16j . The dependence of the tail on r in the in- 
set of Fig. n is consistent with this functional form with 
d = 3 and c ~ 1.3. 



III. FEEDBACK CONTROL OF A BEC 
DROPLET 

The 3D stationary solution of Eq. JHJ is dynamically 
unstable against collapse and expansion. The aim of 
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the present paper is to show that we can stabilize it by 
controlling the scattering length a. When the system is 
about to collapse, we can increase a to prevent the col- 
lapse. Similarly, a decrease in a can prevent expansion. 
Thus, by measuring the density profile of the condensate 
in a nondestructive manner, we can achieve feedback con- 
trol of a to stabilize a BEC droplet. 

An observable quantity in nondestructive phase- 
contrast imaging |13| | is the column density of the con- 
densate, given by 



(a) 0.1 



d co \(x,y,t) = / dz\ip(r,t)\ 



(9) 



where the line of sight is assumed to be in the z direction. 
We use the peak value of the column density, 



D = 4oi(0,0,*), 



(10) 



for the feedback control. 

Let Dq be the target value of the stabilized col- 
umn density. The feedback loop should operate on the 
strength of the interaction g in such a manner that any 
deviation D = D — Dq from the target value Dq will be 
suppressed. We assume that the time derivative of g de- 
pends on D up to the second derivative with respect to 
time: 




g = A(D — Dq) + BD + CD, 



(11) 



where A, B, and C are dimensionless constants. Recov- 
ering the dimensions of D and t by multiplying by N / I 2 
and m£ 2 /h, respectively, where i is an arbitrary unit of 
length, we can write Eq. (fTT|) as 
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47riV 2 m 



A(D ~ Da) + 4^ BD + A^Ti 00 - (12) density profile h/f at t 



FIG. 2: Time evolution of (a) the deviation D = D — Do 
of the peak column density (see Eqs. © and I1UI 1 and (b) 
the fraction of remaining atoms N/No (solid curve) and the 
strength of the interaction g = g — go (dashed curve), for 
A = 10, B = 100, C = 50, At = 0.01, D = 1/tt, and g = 
— (27r) 3 ' 2 . The initial state is the Gaussian wave function ip = 
e~ r / 2 /n 3 ^ 4 , and the value of g starts from (not visible in 
the present resolution). The inset in (b) shows the converged 



100. 



In experiments, a time sequence of phase-contrast im- 
ages is taken by a CCD camera at a certain frame rate 
(e.g., 20 kHz inRef. [3). Therefore, the measurement of 
D is performed at discrete times and the time derivatives 
in Eq. (|llfl are approximately obtained by 



D(t) 

b{t) 



At 
1 

At 2 



[D(t)-D(t-At)}, 



(13a) 



[D(t) - 2D(t - At) + D(t - 2At)], 



(13b) 



where At is the interval between measurements. Accord- 
ingly, g is also changed stepwise as 



g(t) = g(t-At)+At ^A[D(t) - D ] + BD(t) 



■CD(t)}, 

where Eqs. (|13af> and (|13bl) are used for D(t) and D(t). 
The interval At must be made much smaller than the 
characteristic time scale of the dynamics. The At depen- 
dence of the stability is discussed in Sec. El 



We assume that the initial state is the noninteract- 
ing ground state in an isotropic harmonic potential V = 
mw 2 r 2 /2. Henceforth, we take the units of length and 
time to be I = [S/fmw)] 1 ' 2 and uj~ x , respectively. The 
initial state is then described by the Gaussian wave func- 
tion tp — e~ r / 2 /7r 3 / 4 . At t = 0, the trapping potential is 
suddenly switched off and the strength of the interaction 
is set to be g = — (27r) 3 / 2 . The strength of the interaction 
g evolves in time according to Eq. I|14|) . where the value 
of Dq is chosen to be the peak column density of the 
initial wave function, Dq = J dz\ip(x = 0,y = 0, z,t = 
0)| 2 = 1/-7T. The wave function evolves in time according 
to Eq. ©, which is numerically solved by the Crank- 
Nicolson method. We neglect the effect of gravity by 
assuming that it is canceled using, e.g., a technique of 
magnetic levitation [l9| . 

Figure [21 (a) shows the time evolution of D for A = 
10, B = 100, C = 50, and At = 0.01. The value of 
D deviates significantly from Dq only for t < 20, and 
quickly converges to Dq thereafter. This indicates that 



4 



the feedback control of the interaction successfully works 
to stabilize a BEC droplet in 3D. The inset in Fig. [21(b) 
shows the density profile of the condensate at t = 100. 
This converged density profile is found to be related to 
the stationary solution of Eq. JSJ by appropriate scaling. 
The feedback control can thus transform a Gaussian wave 
function into a stationary solution of Eq. Q . During this 
transformation, a certain fraction of atoms are lost from 
the condensate. The solid curve in Fig. [2] (b) shows the 
fraction of atoms remaining within the radius r = 10 
around the center of the condensate, 



N 
iVo 



drAirr 2 ^ 



(15) 



We find that roughly 10% of atoms are scattered away 
from the condensate without returning to the center due 
to the lack of a trapping potential. 

Figure |21 shows the dependence of the dynamics on the 
feedback parameters. Figure|3](a) shows that the param- 
eter A controls the amplitude of the density oscillation. 
The oscillation frequency increases with an increase in 
A, and for large A the system becomes unstable against 
the growth in the amplitude of the density oscillations. 
This is because the term of A in Eq. (| 1 1 ft plays a role of 
pulling the value of D back to Do- If A is too large, this 
pulling causes an overshoot, resulting in oscillations of 
D. Figures 0(b) and (c) indicate that the amplitudes of 
the oscillations decrease with increasing B and C. The 
decay of the oscillations accelerates for larger B [Fig. |3| 
(b)], while it decelerates for larger C [Fig. [31 (c)] . 

Figure 01 shows a stability diagram of the feedback 
control with respect to the feedback parameters B and 
C. We take the value A = 10, since larger values of A 
cause rapid oscillations of the system, as shown in Fig. 01 
(a). We find that the BEC droplet can be stabilized for 
B > 30 and C > -20. 



IV. VARIATIONAL ANALYSIS 

In order to understand in an analytic manner the sta- 
bility of the system subject to the feedback control dis- 
cussed in Sec. IIIII we conduct a variational analysis. 

We employ the Gaussian trial function |2(j 

^o = -J^e-^ + ^\ (16) 



^3/4^3/2 

where R(t) is a variational parameter that characterizes 
the size of the condensate. Equation JSJl is derived by the 
application of the variational principle to the action 



S = / drdtip* 



j— H -\ibr 

dt 2 2 m 



(17) 



Substituting the variational wave function ((16(1 into 
Eq. I |l 7(1 and minimizing the result with respect to R, 
we obtain the equation of motion for R as 



9 



1 



'2tt) 3 / 2 i? 4 ' 



(18) 
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FIG. 3: Dependence of the dynamics of D on the feedback 
parameters A, B, and C. The other conditions are the same 
as in Fig. |2l 



The unstable stationary solution of Eq. 1 (18(1 is obtained 
for R = -g/(27r) 3 / 2 ; hence, if g = -(2tt) 3 / 2 , we have R = 
1. In Fig.Q we plot the Gaussian function ((16(1 with R = 
1 as a dashed curve. Comparing it with the numerically 
obtained stationary solution of Eq. © with g = — (2tt) 3 ^ 2 
(solid curve) , we see that the central density is larger and 
the tail is longer for the stationary solution than for the 
Gaussian function. 

Since R cannot directly be measured in experiments, 
we rewrite the equation of motion in terms of the central 
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FIG. 4: Stability diagram of the BEC droplet under feedback 
control for A = 10. The other conditions are the same as in 
Fig. |21 The open circles show stable points and the filled 
circles show unstable points. 
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FIG. 5: Time evolution of the peak column density D with 
measurement errors given in Eq. (1241 for e — (solid line), 
e = 0.2 (dashed line), e = 0.6 (dotted line), and e = 0.7 (dot- 
dashed line). The feedback parameters are A = 10, B = 100, 
and C = 0. The other conditions are the same as in Fig. [5] 



column density 



D = 2 



dr\ipQ 



ttR 2 ' 



(19) 



We consider small deviations from go = — (2-7r) 3 / 2 and 
Rq = 1, and decompose D and g as D = Do + D and 
g = g + g, where Do = 1/tt. Equation ljl%|l can then be 
rewritten as 

^ = 6 -^' < 2 »> 

where we have kept only the terms linear in \D\ and \g\. 
From Eq. Ijllfl. we have 



g = AD + BD + CD. 



(21) 



Differentiating Eq. I|2U|) with respect to t and using 
Eq. J2U, we obtain 
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If the 3x3 matrix in Eq. (|22|l has an eigenvalue whose 
real part is positive, D diverges exponentially in time. 
Therefore, the condition for stability is that all the eigen- 
values have negative real parts, which can be examined 
by using the Routh-Hurwitz criterion |2lj ] without actu- 
ally solving the eigenvalue equation. The stability con- 
dition of Eq. IL'l'l) is found to be (see Appendix ^] for 
derivation) 



A > 0, 

B > \/27r 5/2 , 

C > 0, 



A < C 



B 



- 1 



(23a) 
(23b) 
(23c) 

(23d) 



These inequalities qualitatively agree with the stability 
diagram in Fig. 0] 



V. EFFECTS OF EXPERIMENTAL 
IMPERFECTIONS ON THE STABILITY OF 
FEEDBACK CONTROL 

We have so far assumed that the peak column density 
D can be measured precisely. However, in real exper- 
iments, there are always imperfections in the measure- 
ment, such as errors and the finite resolution of time and 
space. In this section, we investigate the effects of such 
imperfections on the stability of our feedback control. 

We first study the effects of measurement errors of D 
on the stability. We assume that the measured value with 
an error is given by 



D £ =D(1+ ev rnd ) 



(24) 



where D is defined by Eq. (|1U|> . e describes the error 
level, and v rn d is a random variable that simulates the 
measurement error and is assumed to obey the normal 
distribution e~ 1 ' r " d / 2 /v / 27r. Figure [S] shows some exam- 
ples of the time evolution of D, in which D e is used in 
the feedback equation l|ll|) instead of D. We find that 
the system is tolerant against an error level up to about 
60% in every measurement, and the stability is excellent. 
In Fig. |SJ we see that D has a tendency to decrease with 
an increase in s. This phenomenon is similar to that in 
the case of oscillating interactions 0, Q , where the peak 
density is suppressed by the oscillation of the interaction. 
In the present case, the random fluctuations in g play the 
role of oscillating interactions. 

We next study the effect on the stability of the spatial 
resolution in the measurement of D. We assume that 
due to the finite spatial resolution, the measured value is 
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FIG. 6: Time evolution of the peak column density D with 
spatial resolution given in Eq. 1251 for a = (solid line), 
a = 0.2 (dashed line), a — 0.5 (dotted line), and a — 0.6 (dot- 
dashed line). The feedback parameters are A = 10, B = 100, 
and C = 0. The other conditions are the same as in Fig. [5] 



We have considered the feedback from the peak column 
density D and its time derivatives D and D (Eq. iJHJ), 
and have examined the stability of the system for var- 
ious values of the parameters. We found that stability 
is obtained for a wide range of parameters, as shown in 
Figs. and H 

We have also investigated the stability against ex- 
perimental imperfections, such as measurement errors 
(Fig. [SJ, finite spatial resolution (Fig. |SJ , and finite time 
intervals between measurements. We have found that the 
stability is robust against these imperfections. 

In this paper, we have considered only the simplest 
form of negative feedback. More robust stability may be 
obtained and the stationary state may be reached more 
quickly if more sophisticated methods are used, such as 
the Kalman filter |22j and robust control 23]. 
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filtered by a Gaussian function, 
f 1 

D a = I dxdy-^—^e ^ d co i(x,y), (25) 

where d co \(x,y) is given in Eq. @ and a characterizes 
the spatial resolution. Figure [5] shows the time evolution 
of D, in which D a is used in Eq. Qlljl instead of D. It 
is remarkable that the stability is very robust against 
low resolution. In fact, Fig. HJ shows that the acceptable 
resolution can be almost the size of the condensate itself. 

We have assumed so far that the successive measure- 
ments are performed at an interval of At = 0.01. We have 
also examined the stability for larger At with A = 10, 
B = 100, and C = (with the other conditions the same 
as in Fig. [2J , and found that stability is achieved for a 
time interval of up to At = 0.4. 

As an example, let us consider the case of 85 Rb atoms 
and take the units of length and time to be 3.5/mi and 
16 ms, which corresponds to ui — 10 x 2tt Hz. Then, the 
resolution of the phase-contrast imaging must be a ^ 
1.7/im, and the interval of the measurements must be 
At < 6.4 ms, which corresponds to a frame rate > 160 
Hz. If we use a larger condensate (i.e., a smaller uS), these 
restrictions can be relaxed. 



VI. CONCLUSIONS 

We have shown that a BEC droplet (self-trapped con- 
densate) can be stabilized in 3D free space by the feed- 
back control of the strength of the interaction between 
atoms. By negative feedback on the strength of the in- 
teraction through nondestructive monitoring of the peak 
column density of the condensate, we can prevent the 
condensate from collapsing and expanding. Even start- 
ing from a Gaussian wave function, we can reach the 
stationary state of the GP equation by feedback control. 
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APPENDIX A: STABILITY CONDITION WITH 
THE ROUTH-HURWITZ CRITERION 

According to the Routh-Hurwitz criterion [2l| , all so- 
lutions of a polynomial equation 

a \ n + aiA n_1 H h a n _iA + a n = (Al) 

have negative real parts if (1) all the coefficients a, for i = 
0, 1, • • • , n are real and positive and (2) the determinants 



«1 


a* 


a 5 • 


• 0>2i-l 


a 




CI4 • 


■ 0-2i~2 





<Zi 


a 3 • 


■ fl2i-3 





«0 


a 2 • 











a» 



for i = 2, 3, • • • , n are positive. 

The system described by Eq. I|22|l is stable if all the 
eigenvalues of the 3x3 matrix on the right-hand side have 
negative real parts. The eigenvalue equation is given by 

A 3 + ^-^A 2 +f^— = o. (A3) 

Applying the Routh-Hurwitz criterion to Eq. (|A3|I . we 
obtain the condition for the stability given in Eq. I|23|) . 
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